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C ~ ) , Abstract. Sum of powers l p + ■ • ■ + n p , with n,p g N and n > 1, can be 

expressed as a polynomial function of n of degree p + 1. Such representations 
i are often called Faulhaber formulae. A simple recursive algorithm for com- 

puting coefficients of Faulhaber formulae is presented. The correctness of the 
algorithm is proved by giving a recurrence relation on Faulhaber formulae. 
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1. Introduction 

O . Define f p (ri) = l p + ■ ■ ■ + n p , for p 6 N, n G N + . One can express f p (n) as a 

polynomial function of n of degree p+1. For example, l 2 + • • -+n 2 = ^n 3 + -kn 2 + hn. 
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Such representations are often called Faulhaber formulae, after Johann Faulhaber 
OO ' (1580-1635). In this paper, we study the following problem: 

in 

CO ■ For p £ N, find the Faulhaber formula that represents l p - 

, Let us write f p (n) = a p+1 n p+1 + a n p n + • • • + a\n + a ; it can be proved that 

ao = for all pgN. Clearly finding the Faulhaber formula, given any p6N, can 

ON . 

, be reduced to finding the corresponding coefficients: ax, ■ ■ ■ , a p +i. 

A well-known relation between Faulhaber formulae and Bernoulli numbers, that 
is f p (n) = ^j-j- Xa=o ( P T X ) bin p+1_l with being the i th Bernoulli number (when bi = 
+5), can be used for computing 01, • • • , a p +\. This approach however requires com- 
puting Bernoulli numbers bo, • • • , b p . 

There are various algorithms in the literature for computing Bernoulli numbers. 
These algorithms are generally based on recurrence relations, where bi is computed 
using bo,-- - , bj_i, e.g. see [KB67, AD09]. In this paper, we give a recurrence 
relation on Faulhaber formulae, which yields a direct algorithm for computing the 
coefficients ax, ■ • • , a P +i- 

Structure of paper. In section 2 we give a recursive algorithm for computing coef- 
ficients of Faulhaber formulae. Correctness of the algorithm is proved in section 3. 
Time complexity of the algorithm is also analysed in section 3. 

2. Direct algorithm 

We are interested in computing the coefficients of the Faulhaber formula that 
describes f p (n), for a given p E N. Write f p (n) = a p+ \n p+1 + ■ ■ • + axil. Let us 
consider a table in which rows refer to different values of p, and columns refer to 



the powers of n. The element at the intersection of row i and column j, denoted 
a(ij), is meant to represent the coefficient of w? in the polynomial describing fi(n). 
Sec figure 1. Note that elements at with j > i + 1 arc all zero. 
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Figure 1 . Intuitive description of the direct algorithm 

Our goal is therefore to find the numbers in the row corresponding to p. We 
proceed inductively: First, row is filled, then we fill row 1, . . . , till the row 
numbered with p is filled. The algorithm starts with placing 1 at position (0, 1) in 
the matrix, thus filling row 0. This step reflects fo(n) = 1° + • — h n° = n. In order 
to fill the i th row, with i > 0, we follow the rules below: 

(1) For 1< j < i + 1, let a (lJ) = ja^j^y 

(2) Next, we compute a^i) as a^i) = 1 — Xa<j<i+i a (i-j) ■ differently, 
a(i.i) is chosen such that the sum of the numbers that appear in each row 
equals I. 

The procedure stops when the row corresponding to p is filled. Below, it is proved 
that a(p j) , for 1 < j < p+ 1 , is the coefficient of v? in the polynomial that represents 
f p (n). Algorithm 1 implements this procedure. 

Note that, since filling row i only requires elements of row i — I, the algorithm 
only stores a vector, instead of the matrix of figure 1. 

Algorithm 1 Computes coefficients of Faulhaber formulae 
Require: peN 

ai := I 

for (i := 1; i < p; do 

s:=0 

for (j := i + 1; j > I; j ) do 

s := s + <ij 
end for 
ai := 1 — s 
end for 

return a\, ■ ■ ■ , a p+i 

Using the table of figure 1 and the presented algorithm, we make the following 
simple observations about coefficients of Faulhaber formulae. 
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• The coefficient of n p+1 in f p (n) is for any p E N. This can be proved 
by induction: a (PiP+ i) = ^-a (p _i iP ) and a ( i >2 ) = |. 

• The coefficient of n p in / p (n) is |, for any p e N. This can be proved by 
induction: a (PiP ) = |a (p _i iP _i) and a (M) = |. 

• The coefficient of n p ~ 2 in / p (n) is zero, for any p > 3. This can be proved 
by induction: a( PiP _ 2 ) = p^2 a ( P -i,p-3) and 0(3,1) = °- 

3. Recurrence relation on Faulhaber formulae 

In this section we prove that algorithm 1 correctly computes the coefficients of 
Faulhaber formulae. For this, first, a recurrence relation on Faulhaber formulae is 
proved. 

Lemma 1 (Recurrence on Faulhaber formulae). 

• fo(n) = n, 

• forp> 0, f p (n) =pfa f P -i(t) dt + (1 - p J* f p -i(t) dt)n 

Proof. The first part, fo{n) = 1° H + 2° = n, can be proved by straightforward 

induction. In the following, therefore, we assume p > 0. To prove the second part, 
we recall the following relations (e.g. see [AS72, chapter 23]): 

I fp-i(n) = \{B p (n + 1) - B p (0)), for p g N + . 
II S h a Bi{t)dA = j^(B t+1 (b) B i+1 (a)). 
Ill Bi(n+1) -Bi(n) = in 1 - 1 , for n E N,i > 1. 

where Bi{t) is the i th Bernoulli polynomial. Using (I) and (II), we get: 

/>/ P -i(i) dt 

= f a B p (t + l)-B p (0)dt 

UBp+^b + 1) - B p+1 (a + 1)) - B p (0)(b - a) 



P+i 



Therefore 



l-foPf P -i(t)dt 

= l--l T (B p+1 (2)-B p+1 (l))+B p (0) 
= B p (0) 

The last simplification step is due to (III). Note that since p > 0, we have p+1 > 1, 
satisfying the precondition of (III). Similarly 

foPf P -i(t) dt 

= -^{Bp+iin + 1) - B p+1 (l)) - nBp(0) 



As a result 



IoPfp-i(t) dt+(l- /oP/ P -i(t) dt)n 
= ^ T (B p+1 (n+l)-S p+1 (l)) 
= F i T (B p+1 („+ 1)-S P+1 (0)) 



The last simplification step is again due to (III): B p+1 (l) — B p+ i(0) = 0, with 
p + 1 > 1. Finally 

= -L I (B p+1 (n + l)-B p+1 (0)) 
= JoPfp-iit) dt+(l- JoPfp-i(t) dt)n 
This completes the proof. □ 

Now, we are ready to prove the correctness of algorithm 1. 

Theorem 1 (Correctness). Given p 6 N. algorithm 1 outputs the coefficients of 
the Faulhaber formula that represents f p (n). 

Proof. Let us assume the coefficient of n? is a in fi(n), for some 1 < j < i + 1, and 
the coefficient of n^ 1 in /j_i(n) is (3. From the recurrence relation of lemma 1, we 
get a = 4/3. This is simply because J™ t k alt = -j^p[n k+1 , where k is a any positive 
rational number. This directly results in the way algorithm 1 recursively computes 
a(ij), for 1 < j < i + 1, and i > 0. 

Now, note that f p (l) = 1, for any p e N. Moreover, note that / p (l) = a p+ \ + 

h a\. Therefore, a\ = 1 — (a p+ i + • • • + a 2 ). This immediately results in the way 

algorithm 1 computes !), for i > 0. □ 

Below, we turn to time complexity of algorithm 1. To measure the computational 
complexity, we count the number of multiplication and addition (or, subtraction) 
operations that are performed on rational numbers. Assignments to constants, and 
incrementing and decrementing natural numbers (i.e. counters in the algorithm) 
are thus assumed to take negligible time. 

Theorem 2 (Time complexity). Time complexity of the direct algorithm is qua- 
dratic in p. 

Proof. Note that the outer for loop is repeated p times, and the inner for loop is 
repeated 1 + ■ ■ ■ + p = \p{p + 1) times. It is straightforward to see that, for p 6 N, 
the number of addition operations is \p{p+ 1) +p, and the number of multiplication 
operations that are performed on rational numbers is \p(j> + 1). □ 
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